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We study the simplest irreversible ballistically-controlled reaction, whereby particles having an 
initial continuous velocity distribution annihilate upon colliding. In the framework of the Boltzmann 
equation, expressions for the exponents characterizing the density and typical velocity decay are 
explicitly worked out in arbitrary dimension. These predictions are in excellent agreement with the 
complementary results of extensive Monte Carlo and Molecular Dynamics simulations. We finally 
discuss the definition of universality classes indexed by a continuous parameter for this far from 
equilibrium dynamics with no conservation laws. 



Systems with reacting particles model a rich variety of 
phenomena and provide prominent situations to develop 
and test the foundations of non-equilibrium statistical 
mechanics. In this context, the diffusion controlled first 
order annihilation process (A + A — > 0) has been ex- 
tensively studied and the corresponding decay kinetics is 
well understood. On the other hand much less is known 
in the contrasting case where the reactants move ballis- 
tically between the collision events, despite the relevance 
of such motion for growth and coarsening processes ^) . 
A few theoretical results are available in d = 1 dimension 
for such irreversible kinetic processes with discrete initial 
velocity distributions. In a pioneering work, Elskens and 
Frisch show from combinatorial considerations that the 
particle density n{t) decays like \j\fi for the simplest bi- 
nary velocity distribution || . Powerful generalizations of 
this result were obtained still in ID, either for a larger 
class of stochastic ballistic annihilation and coalescence 
models ^) or from kinetic theory for discrete multi- 
velocity distributions || FS. No exact results could be 
obtained for the generic case of continuous distributions, 
where the decay exponents have been computed numer- 
ically H |], Recently however, Krapivsky and Sire 
considered the latter situation in the framework of the 
Boltzmann equation (relying on the so-called "molecular 
chaos" factorization O]) and derived bounds for the ex- 
ponents as well as their leading large d behavior. The 
existing body of literature has essentially focussed nu- 
merically on the one dimensional case, and no accurate 
predictions seem to be available for the decay exponents. 

In this Letter, we obtain predictions for the decay ex- 
ponents and velocity distribution (assumed initially con- 
tinuous) , revisiting Boltzmann kinetic theory in arbitrary 
dimension, with the explicit inclusion of non Gaussian 
corrections to velocity distributions. These predictions 
are compared both with the existing numerical results 
in ID and the expressions derived in |8[ |o), and further 
tested against extensive numerical simulations in dimen- 
sion 2 and 3, following two complementary routes: we 
first solve the mean-field non-linear Boltzmann equation 
describing the annihilation process by means of a Monte 
Carlo scheme, which validates the analytical expressions 



obtained within the molecular chaos framework; second, 
we go beyond mean-field and investigate the exact de- 
cay kinetics by implementing Molecular Dynamics simu- 
lations. The two numerical approaches yield the same ex- 
ponents, in excellent agreement with the analytical pre- 
diction. Finally, we address the question of universality 
in this process || by partitioning the possible contin- 
uous velocity distributions into groups associated with 
the same asymptotic dynamic scaling behaviour, akin to 
equilibrium universality classes. 

We consider an assembly of identical spherical particles 
with radius a in dimension d, with initial velocity distri- 
bution f(v, t = 0) and random initial positions. Particles 
follow free flight motion until a collision occurs which re- 
sults in the removal of both partners. We are interested 
in the time evolution of density n(t) = j f(v,t)dv and 
typical velocity v(t), related to the kinetic temperature 
T(t) defined as the variance of the velocity distribution 



T(t) = -J-r 

V ' n(t) 



v 2 f(v,t)dv = (tJ) 2 



(1) 



Insight into the decay kinetics may be gained by writing 
the rate equations for n and T 
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= —Lu(t)nT co u — —au)(t) nT, 
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where the first line stands for a definition of the instanta- 
neous mean collision frequency uj, while r co ii is the time 
dependent total kinetic energy of a colliding pair, which 
is thus dissipated in a binary encounter, as stated by 
the rhs equality in Eq. (^). On dimensional grounds, 
the collision frequency is expected to scale like the in- 
verse time, which together with Eqs. (||) and (J3j> implies 
an algebraic time decay for n and v, as well as a time- 
independent energy dissipation parameter a [defined in 
Eq. (||) as a — T co \\/T]. We therefore introduce two 
exponents £ and 7 such that n(t) oc t~^ and v oc 
(and T oc i -27 ). With a ballistic dynamics controlled 
by the mean- free-path I oc l/(n<7 rf ~ 1 ), the collision fre- 
quency may be written as the ratio v/i. From uj oc 1/t, 
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we obtain the scaling relation £ + 7 = 1 0, ^, 
which may be combined with the ratio of Eqs. 
(§) to give 
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and 
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7 
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(4) 



Since particles with a higher velocity are likely to disap- 
pear with a higher rate than the average particle with 
temperature T, we expect a = T co \\/T to be larger than 
1, so that the typical velocity should decrease with time 
[7 > from Eq. (Q)]. This moreover explains the failure 
of the naive mean-field picture where the density decay 
rate is written n oc — n 2 , so that n(t) oc 1/t. This trans- 
parency limit would hold in the absence of collisional cor- 
relations (a = 1) which becomes only asymptotically ex- 
act in the limit of infinite dimension d. 

We now turn to the computation of a within the molec- 
ular chaos framework, which is a priori an uncontrolled 
approximation. It will however be shown to capture 
the essential collisional correlations missed by the naive 
mean- field argument, and to provide decay exponents in 
excellent agreement with their numerical counterparts. 
The corresponding Boltzmann equation reads 



df(v,t) 
dt 



= -f(v,t) J dw\v-w\f(w,t), 



(5) 



which implies that if the initial distribution behaves like 
a power law near the velocity origin, this property 
is preserved at subsequent times by the dynamics, which 
in turn should affect the exponents £ and 7, expected 
to depend explicitly on /1 (as appears on the analytical 
predictions of Ben-Nairn et al J§ £ = (2d+2(i)/(2d+2u + 
1), or on the bounds derived by Krapivsky and Sire [10|). 
Looking for a scaling solution of the kinetic equation^) , 
we introduce a rescaled velocity c = v/v and rescaled 
single particle distribution function tp through 



f(v,t) = ^ v>(c,t), 
v 



(6) 



so that tp(c, t) is the probability distribution function 
of the velocity c at time t, satisfying the constraints 
J (pdc = 1 and J c 2 tpdc = 1 at any time. If f(v,t) 
evolves into a self-similar decay state, the only relevant 
time dependence occurs via n(t) and v(t), so that tp(c,t) 
no longer depends on time and the evolution equation for 
/ (assumed isotropic) translates into 
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C12 



(7) 



where ((...)) = /(. . .)tp(ci)tp(c 2 )dcidc2 so that (ci 2 ) = 
(\c 1 — c 2 \) denotes the rescaled collision frequency. 

Equation (j^) may be considered as an eigenvalue prob- 
lem for a, which has been computed numerically in ID 



p3 . However, it is useful to reformulate Eq. (fa) into 
an infinite hierarchy of consistency relations obtained by 
computing the corresponding moment of order p: 
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Note that the special case p = 2 coincides with the defi- 
nition of a through the kinetic energy dissipation as ex- 
pressed by Eq. §): a - T col] /T = (c 12 c?)/((c 12 )(cf)). 
We look for explicit solutions by expanding tp in a basis 
of Soninc functions Jl3l 



ip(c) = M{c) 



E 



a n Sn( c ) 



(9) 



where the polynomials S n are orthogonal with respect to 
the Gaussian weight M. (c) . Computing the averages in- 
volved in from the functional expression (Q) provides 
a system of equations for the coefficients a n . 

In practice, only a few terms are required in the ex- 
pansion (^|) in order to get a precise estimation for a, 
provided relations of lowest order p as possible are re- 
tained among the hierarchy (||). In this respect, taking 
the limit of vanishing velocity of (0) yields the "optimal" 
relation involving a and moments of tp of order 1: 



1 - 



(C12) 



(10) 



that we consider as the first equation of (g) correspond- 
ing to the limit p — > + . At Gaussian order for tp [i.e. 
truncating (|^) at order n = 0], it is straightforward to 
get 
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1 - 



V2 



(11) 



which, together with Eq. 
timation for £: 

£0 



yields the zeroth order es- 



2d+2fi 



2(d + A* + l) - V2' 



(12) 



It is noteworthy that in the limit of large dimension, we 
obtain £ ~ 1 - cT^l - I/V2) + 0(l/d 2 ) irrespective of 
/x, which has been shown to be the exact l/d behaviour 
within Boltzmann molecular chaos framework fl(i[| . The 
first non Gaussian correction is carried by a 2 {a\ iden- 
tically vanishes from the definition of temperature | fl4| ) 
and this coefficient is related to the kurtosis of the veloc- 
ity distribution: a 2 is proportional to the fourth cumu- 
lant (cf) — 3(cf ) 2 , where Cj is a given Cartesian coordinate 
of c. After a lengthy calculation performed at linear or- 
der in a 2 , we obtain 
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a 2 



d(3 - 2V2) 



4d 2 + 6/2 + d(6 + 4// - V2) 
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Q0 + T6d fl2 - 



(13) 



(14) 
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The above predictions rely on a perturbative expan- 
sion starting from the Maxwellian M. (regular at v = 0) 
and are therefore expected to be particularly relevant for 
fi close to 0. The agreement with the existing numerical 
data is excellent; an accurate estimation has been re- 
ported in ID within molecular chaos for the much stud- 
ied = case [10 : £ = 0.769(5) whereas we obtain 
at zeroth order £o = 0.773 from ([l2]) and at second or- 
der £ 2 = 2/(1 + a 2 ) = 0.769(3) from Eq. @. This 
exponent is compatible with its counterpart extracted 
from the exact dynamics (0.78 in M). Moreover, we 
have investigated numerically the annihilation dynamics 
in higher dimensions by means of a) the Direct Simu- 
lation Monte Carlo procedure jl5) (DSMC) solving the 
non-linear homogeneous Boltzmann equation M) and b) 
Molecular Dynamics simulations (MD) implementing the 
exact dynamics with periodic boundary conditions [fL6[ . 
The DSMC technique provides precise data for the ve- 
locity distributions and decay exponents, and allows to 
test the validity of the analytical truncated expansion of 
the scaling form ip, leading to ( jll| ) or (|l~J). Alternatively, 
MD results assess the reliability of the molecular chaos 
ansatz, but are more demanding on computer resources: 
on the one hand, the system needs to reach very low 
densities in order to develop the self-similar decay stage 
where f(v,t) takes the scaling form (^), but on the other 
hand, the mean free path I which increases with time like 
IP must remain smaller than the simulation box size L, 
which provides a lower bound for n(t) or equivalently an 
upper bound for accessible times before finite size effects 
hinder the precise determination of £ and 7. In practice, 
we considered systems with TV = 10 5 — 5 10 5 particles in 
MD and N = 10 6 — 10 7 in DSMC where it is further 
possible to average over 10 3 to 10 4 replicas to increase 
the statistics of the velocity distributions, which is cru- 
cial for computations at large times with a concomitant 
low number of particles left. 



(ID) values of fi 


-4/5 -1/2 


Prediction [H, E(| 

Simulation [3] 
f 2 from Eq. ( L44) 


0.28 0.5 0.666 
0.32/0.37 0.56/0.60 0.769 
0.32 0.60 0.769 


TABLE I: Decay exponent £ in 1 dimension. 


(2D) values of [i 


-1 -1/2 3 


Prediction [§, 

Simulation 
§2 from Eq. (p|) 


0.66 0.75 0.800 0.91 
0.75 0.83 0.870 0.97 
0.76 0.84 0.870 0.95 



TABLE II: Exponent £ in 2D; the simulation data are the 
Monte Carlo results of the present work. 

The results of two dimensional simulations are shown 
in Fig. |Tj where it appears that the MD data are 
fully compatible with DSMC, although less precise. For 
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FIG. 1: Evolution of the density (lower sets) and kinetic tem- 
perature (upper sets), normalized by their initial values. At 
t — 0, the velocity distribution is Maxwellian (fi = 0), with 
a collision frequency denoted uio. MD results are shown by 
symbols (circles for n and crosses for T) and DSMC by con- 
tinuous curves. The dashed lines have slopes given by the 
theoretical predictions. Inset: check of the scaling relation 
£ + 7 = 1 where ny/T is expected to scale like t~ ?-7 ; the 
dashed line has slope -1. 



uiot ~ 10 3 , the MD density and temperature tend to sat- 
urate, which corresponds to the upper time limit where 
£ ~ L, and the subsequent evolution is discarded. The 
predictions £0 = 0.872 and £2 = 0.870 for n — (indistin- 
guishable in Fig [l]) are in good agreement with the sim- 
ulations, irrespective of the initial f(v) chosen (we con- 
sidered several distributions with the constraint /i = 0, 
see below the discussion concerning universality). The 
above exponent is compatible with that reported in the 
context of a multi-particle lattice gas method (0.87 p7[). 
Moreover, the initial spatial configuration is irrelevant 
(the long time dynamics and rescaled velocity distribu- 
tions are the same starting from a fluid-like structure or 
from various crystalline arrays), and the scaling relation 
£ + 7 = 1 is seen to be well obeyed in the asymptotic 
regime (inset of Fig. |]). The same scenario holds in di- 
mensions 3 and 4, where the predictions at zeroth and 
second order are very close, and indistinguishable from 
the numerics (£0 — 0.91 in 3D and 0.93 in 4D for /j, = 0). 
However, the agreement is expected to become worse as 
fi deviates from (with /1 > — d to ensure proper nor- 
malization). This is confirmed in Tables | and |H| which 
summarize the results obtained for various /1, with com- 
parison to the theoretical prediction of Ben-Nairn et al. 
[ 8| (coinciding with the lower bound for £ obtained in 
[HI, the upper bound being 1). For [i = 0, the non 
Gaussian parameter a 2 is small [with an even smaller 
correction to a due to the prefactor v / 2/(16cf) in (p4^)]. 
This fourth cumulant however rapidly increases with /x, 
so that inclusion of higher order terms [n — 3. . . in (^)] 
would be required to obtain the same level of accuracy 
as for regular distributions near the velocity origin. 
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FIG. 2: Plots of ip(ci) I M(ci) versus a in 2D. The inset shows 
4 different initial distributions with [i — 0, one of them be- 
ing Gaussian [thus corresponding to the flat curve (circles)]. 
These distributions having very different a 2 at t = collapse 
onto a master curve in the asymptotic scaling regime (main 
graph). The thick curve is the prediction I + a 2 S2{c 2 ) where 
a 2 is given by Eq. (|f|) and S 2 {x) = x 2 /2 ~ 3x/2 + 3/8. The 
symbols (stars, crosses, pluses and circles) refer to the same 
distributions at late times (main graph) and at t = (inset). 
The results have been obtained by averaging over 10 4 replicas 
of a system with N = 5 1 6 particles. 

In the remainder, we consider the possibility to define 
universality classes for ballistic annihilation kinetics, in 
the following sense: does /i completely specify the asymp- 
totic velocity distribution and decay exponents, irrespec- 
tive of further details concerning the initial conditions 
P|? To answer this question we have run several simula- 
tions (MD and Monte Carlo) corresponding to different 
initial conditions sharing the same \x, for several values 
of this parameter. The corresponding decay exponents £ 
and 7 are monitored, which provides a first test, however 
quite insensitive to possible non Gaussianities (see above 
the numerical proximity between £o an d the non Gaus- 
sian corrected £2)- A more sensitive and severe probe is 
provided by the kurtosis et2, which may be computed in 
two different ways: first from its definition involving the 
fourth cumulant (cj) — 3(c|) 2 , or alternatively from the 
direct computation of ip(c)/M(c), which may further be 
compared to the analytical expansion 1 + 026*2 (c 2 ) with 
a 2 given by Eq. ( |l3| ) (recall that ai = 0). The latter 
method is illustrated in Fig. |^ where the four initial 
distributions shown in the inset evolve after a transient 
towards the same attractor, that is furthermore in quan- 
titative agreement with the Sonine prediction. Moreover 
the same values of £ and 7 are measured within statistical 
inaccuracy for the 4 distributions. We have observed the 
same phenomenology for /i 7^ 0, which points to the rele- 
vance of defining universality classes of initial conditions 
as distributions having the same regularity exponent /z, 
as conjectured in ID for ji = Q. 

In conclusion, we have shown that the non trivial 
dynamic scaling behaviour of ballistic annihilation may 
be investigated within Boltzmann kinetic theory, and 
accurate decay exponents have been explicitly worked 



out. Their evaluation (g_2|) at zeroth order turns out 
to be straightforward, but follows from a kinetic equa- 
tion and is therefore specific to the precise model consid- 
ered here. A more versatile approach that would apply 
to any ballistically controlled reaction (including coales- 
cence with arbitrary conservation laws, with or without 
stochasticity in the reactions) consists in reconsidering 
the rate equations (J2j> and ||), and identify the proper 
energy dissipation parameter a before approximating it 
assuming a Gaussian velocity distribution. This "model- 
independent" approach gives a = 1 + l/(2d) in the par- 
ticular case of pure annihilation, which corresponds to 
£ = 4d/(4d + 1) (i.e. 0.8, 0.89 and 0.92 in dimensions 
1, 2 and 3) in reasonable agreement with the exponents 
mentioned above (0.77, 0.87 and 0.91 respectively). We 
conjecture that the exponent £ = 4c£/(4d+l) becomes ex- 
act when the particles annihilate with probability p (and 
collide elastically otherwise) , in the limiting case p — > + 
(whereas p = 1 for "pure" annihilation). This hopefully 
provides an illustration of the central role played by the 
energy dissipation parameter a in ballistically controlled 
reactions, and calls for further investigations with more 
involved reactions. 
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